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I. INTRODUCTION 



Many approaches in physical chemistry deal with model systems consisting of clas- 
sical rigid molecules. In studying such systems by the method of molecular dynamics 
(MD), three problems arise at least: (a) the choice of suitable parameters for describing 
a state of the system in phase space, (b) the application of an efficient algorithm to inte- 
grate numerically the equations of motion, and (c) the exact conservation of the rigidity 
of molecules during an approximate integration. The question of how best to handle these 
problems is one which has been keenly debated and the relative merits of a number of 
various schemes have been devised. 

The molecular approach treats dynamics of the system in view of translational and 
rotational motions. In the classical scheme [1] , three Eulerian angles are used to represent 
the same number of rotational degrees of freedom of the molecule. A numerical integra- 
tion of the corresponding equations of motion was performed in early investigations [2, 3]. 
It has been soon established [4, 5] that this integration is very inefficient because of sin- 
gularities whenever the azimuthal angle of the molecule takes a value of or ir. Although 
the singularities can be avoided by applying different sets of Eulerian angles, this requires 
complex manipulations with time-consuming trigonometric functions. In singularity free 
schemes, the orientations of molecules are expressed in terms of either quaternions [6-10] 
or principal-axis vectors [6]. The last scheme has been derived extending the symmetry 
vector method [11, 12] for diatomics to an arbitrary rigid body. 

In the atomic approach [13], the phase trajectories are considered as translational 
displacements of individual molecular sites. Such particles move independently under the 
potential-energy forces and constraint forces, introduced to hold inter-atomic distances 
constant. This approach is intensively exploited in MD simulations since usual algorithms 
for integration of translational motion can be applied here. However, the atomic technique 
is sophisticated to implement for point molecules and when there are more than two, three 
or four interaction sites in the cases of linear, planar and three-dimensional molecules, 
respectively, because then the orientations can not be defined uniquely [14]. Moreover, to 
reproduce exactly the rigid molecular structure for arbitrary polyatomics, it is necessary 
to solve complicated systems of six nonlinear equations per molecule at each time step of 
the integration process. 
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Usually, high-order predictor-corrector algorithms [15, 16] are applied to integrate the 
equations of rotational motion [2, 7, 8]. Such algorithms, being very accurate at small 
time steps, quickly become unstable and can not be used for greater step sizes [14]. Small 
time steps are impractical in calculations, because too much expensive computer time 
is required to cover the sufficient phase space. At the same time, translational motion 
is successfully integrated with the lower-order Verlet [17], leapfrog [18], velocity Verlet 
[19] and Beeman [20] algorithms, owing their simplicity and exceptional numerical stabil- 
ity (for example, the equations of atomic motion are integrated within the usual Verlet 
framework [13, 14]). However, original versions of these algorithms were constructed on 
an assumption that acceleration is velocity- independent, and, therefore, they can not be 
applied directly to rotational dynamics. Analogous pattern arises with translational mo- 
tion in the presence of external magnetic fields or when relativistic effects are important 
and it is necessary to take into account internal fields of moving charges. 

To remedy that omission, Fincham [21, 22] has proposed explicit and implicit versions 
of the leapfrog algorithm for rotational motion in which angular momenta are involved 
into the integration. In the case of a more stable implicit version, this leads to a system 
of four nonlinear equations per molecule for the same number of quaternion components, 
which is solved by iteration [22]. Ahlrichs and Brode have derived a hybrid method [23] 
in which the principal axes are considered as pseudo particles and constraint forces are 
introduced to maintain their orthonormality. The evolution of principal axes in time can 
be determined using a recursive solution for exponential propagators. In such a way some 
difficulties of the cumbersome atomic technique have been obviated. But the algorithm 
is within the Verlet framework and does not involve angular velocities. Therefore, it 
is impossible to extend it to a thermostat version or to an integration in the presence 
of magnetic fields. Moreover, the pseudo-particle formalism does not contain molecular 
torques, so that it is not so simple matter to apply it to systems with point multipoles. 
Finally, the recursive method [23] as well as the rotational-motion leapfrog algorithms 
[22] appear to be much less efficient with respect to the total energy conservation than 
the atomic-constraint technique. 

In the present paper we develop the idea of using principal-axis vectors as orientational 
variables. We involve the velocities and molecular torques explicitly and show that the 
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rigidity problem can easily be resolved in our approach without any additional transfor- 
mations. The paper is organized as follows. The equations of motion for orientational 
matrices are obtained in Sec. 2. The question of how to integrate these equations within 
the Verlet framework in velocity form is considered in Sec. 3. A comparison of different 
approaches, based on actual MD simulations of water, is presented in Sec. 4. Concluding 
remarks are added in Sec. 5. 



II. EQUATIONS OF MOLECULAR MOTION 

Let us consider a system composed of N identical rigid molecules with M atoms. 
We split evolution of the system in time t into translational and rotational motions. The 
translational motions are applied with respect to the molecule as a whole and can be 
described by the 3N (i — 1, . . . , N) Newton equations 

J 2 N;M 

m-^ = (1) 

where r*j = J2^f rn a r^/m and rf are the positions of the centre of mass and atom a of 
molecule i, respectively, m = m a and m a denote the masses of a separate molecule 
and partial atoms, and are the atom-atom forces between two different molecules. 

To analyze rotational motions, we introduce the sets e = (ei, 62,63) and u % = 
(u\,U2,u % 3 ) of orthogonal unit vectors characterizing the laboratory fixed coordinate 
system, L, and the moving coordinate system, S\ attached to molecule i, respectively. 
Orientations of the S*-system with respect to the laboratory frame can be defined as 
mJ, = or merely u l+ = Aje+, where e+ and u l+ are vector-columns, 

a l a g = u l a -e.p are components of the rotational matrix Aj and a,f3 — 1,2,3. Let us 
place the origin of the S*-system in the centre of mass of the i-th molecule and direct 
the axes of this system along the principal axes of inertia. The principal components of 
angular velocities, fli = Q\u\ + Q^u^ + f}\u\, obey 3iV Euler equations [2], 

Ja^f = K(t) + (jp - j,)ni(t)n;(t) , (2) 

where (a, (3, 7) = (1, 2, 3); (2, 3, 1) and (3, 1, 2). Here Ji, J 2 and J 3 are the independent on 
time principal moments of inertia of the molecule, J2j-a,b^i'xF"ij — ^i e i + ^2 e 2 + k l 3 e 3 = 

4 



K\u\ + K l 2 u 2 + K\u\ is the torque exerted on molecule % with respect to its centre of mass 
due to the interactions with the other molecules, Kf = A^kf ', where K L = (K{, K 2 , K\), 
ki = (k\, k\, k\) and 81 = r°— ty Let A a = (A", Af , Ag) + be a vector-column of positions 
for atom a within the molecule in the S*-system, i.e., 81 = /S.\u\ + A 2 u 2 + A3W3. By 
construction of the S*-system the conservative set (a = 1, . . . , M) of vectors A a is the 
same for each molecule and defined by its rigid geometry. Then the positions of atoms 
in the L-system at time t are r"(t) = r^t) + A+(t)A a , where A+ denotes the matrix 
transposed to A. 

Usually, the elements of orientational matrices Aj are expressed via three Eulerian 
angles which can be chosen as follows: cos#, = e 3 -w 3 , cosy?i = e 2 -(e 3 xu 3 )/|e 3 X'U3| and 
cosipi = U2«(e 3 Xw 3 )/|e 3 Xit 3 |. Then principal components of angular velocity are Q\ = 
9i sin ipi — (pi sin 6i cos ipi, Q\ = 6i cos ipi + (pi sin 6i sin ipi and i? 3 = pi cos 6i + ipi- As was 
mentioned earlier, the equations of motion are singular in this case. The most notorious 
demonstration is the expression (Q l 2 sin-^j — Q\ cos^)/ sin^j for the generalized velocity 
(pi from which it follows that pi — > 00 when 9i tends to zero or n. This leads to serious 
technical disadvantages for the application of Eulerian angles to numerical calculations. 
It is worth mentioning that the rigidity of molecules is conserved automatically in this 
approach, i.e., |^(t)| 2 = (A+(t)A a ) + (A+(t)A a ) = A a+ A l (t)A+(t) A a = |A a | 2 , where 
the property AA + = I of rotational matrices has been used and I is the unit matrix. In 
other words, the matrix Aj remains an orthonormal one for arbitrary values of Eulerian 
angles. 

As is now well established [6, 7], at least four orientational parameters per molecule 
must be used to avoid the singularities. In this case the matrix Aj is a function of these 
parameters which constitute the so-called quaternion q i = (£j, rji, Q, Xi)- It is necessary 
to emphasize that the matrix Aj is orthonormal if the quaternion satisfies the equality 
Qi — Q + Vi + Ci + Xi — 1- I 11 practice, however, the equations of motion are not solved 
exactly, so that this constraint will only be satisfied approximately. The simplest way to 
achieve the required unit norm at all times lies in multiplying each quaternion component, 
associated with the same molecule, by the common factor 1 j \J~q? at every time step of 
the numerical integration [7, 22]. 

In the mentioned above approaches, orientations of the S*-system with respect to the 
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laboratory frame were defined by the matrix Aj, where each of the nine elements \a l a/3 \ < 1 
of which is a function of either three Eulerian angles or four quaternion components. Now 
involving no Eulerian angles and quaternions, we merely consider all these elements as pa- 
rameters which represent the rotational degrees of freedom. The elements a l a g are, in fact, 
the Cartesian coordinates of principal axes of molecules in the laboratory frame. They 
are not independent as it follows from the requirement AjA+ = I imposed on rotational 
matrices. For example, the first three elements a l n , a\ 2 and a\ 3 can be expressed via the 
rest of others from the vector relation u\ = u l 2 Xu l 3 , reducing the number of orientational 
parameters per molecule from 9 to 6. The remaining six elements are connected by three 
constraints, namely, u\-u\ = 1, u 3 -u % 3 = 1 and u\'U\ = 0. Thus, among these six ele- 
ments we can choose arbitrarily three ones, not belonging the same row of the matrix, 
to form an independent set, but only with a few exceptions. Indeed, let a 21 , a 22 and 033 
be chosen as independent elements and one considers a particular case, when a 33 = ±1. 
Then from the equality u l 3 -u 3 = 1 it immediately follows that a 31 = a 32 = 0. From the 
next equality u l 2 -u 3 = we find that a 23 = and, finally, the third relation u 2 -u l 2 = 1 
yields the constraint a 21 + a 22 =l concerning the variables which were assumed to be as 
independent quantities. 

The reason of this situation is similar to that existing in the case of using Eulerian 
angles, where the singularities have appeared at Qi = or n, i.e., at a 33 = cos#j = ±1. 
It indicates again about the impossibility to derive singularity free equations of motion 
involving only three orientational variables per molecule. As was pointed out earlier, 
four orientational parameters avoid the singularities for arbitrary polyatomic molecules. 
Nevertheless, a larger number of parameters can also be acceptable, but this leads to an 
increased number of constraints. For instance, there is one constraint per molecule for 
quaternion variables, while there are three or even six constraints for six or nine parameters 
in our case. From this point of view, such an original presentation [6] of the matrix 
approach has no advantages with respect to the quaternion method. However, as we shall 
show in the next section using some specific properties of the matrix representation, the 
constraints can be satisfied intrinsically within particular integration frameworks without 
any additional transformations. 

The equations of motion for dynamical variables a l a/3 can be found as follows. From 
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the definition du' l a /dt = f2iXu l a of angular velocity and the orthonormality of sets e and 
u\ we obtain a\ a = n\a* 2a - Q\a\ a , a 2a = Q\o^ a - Q\d la and d| Q = Q\a\ a - Q\a 2a , or in 
the matrix form 



(3) 



where 



a = 



Q\ -Q\ 
-Q\ Q\ 



\ 



a 



i -n\ o J 



(4) 



are antisymmetric matrices associated with angular velocities Then differentiating 
relations (3) with respect to time, one obtains the 9iV (i — 1, . . . , N) scalar equations of 
motion 



(5) 



where fli are defined according to Euler equations (2) and angular velocities are excluded 
from equalities (3), i.e., fij = AjA^~. In such a way we construct the coupled set (1), (5) 
of 12 N differential equations of type ^({r*j, r i: Aj, Aj, Aj}) = in terms of the 12N gen- 
eralized coordinates {r^, Aj}. If an initial state {r i (t ),r i (t ),A i (t ),A i (t )} is specified, 
the time evolution {ri(t), Aj(t)} of the system can be unambiguously determined by (1) 
and (5). 



III. INTEGRATION WITHIN THE VELOCITY VERLET FRAMEWORK 

The equations of motion obtained must be complemented by an integration al- 
gorithm in order to be applicable for actual simulations. As was demonstrated for the 
atomic approach [13, 14], a very efficient technique follows from the Verlet algorithm. 
The same framework has been used in the pseudo-particle formalism [23]. However, the 
Verlet algorithm in its original form [17] does not involve velocity explicitly into the 
integration process and, therefore, it can not be applied to equations of motion with 
velocity-dependent accelerations, as in our case (see eq. (5)). Because of this we shall 
work within a velocity form [19] of the Verlet method. 
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A. Basic ideas 



Let {ri(t ), Viito), Aj(t ), Aj(t )} be a spatially-velocity configuration of the system 
at time to- On the basis of equations (1) for translational motion we can calculate the 
translational accelerations fj(t ) using molecular forces Fi(t ). Then, according to the 
first line of the velocity Verlet integrator, the positions of the centres of mass of molecules 
(i — 1, . . . , N) at time t + At are 

Ti(to + At) = ri (to) + fi(t )At + r l (t )At 2 /2 + 0(At 3 ) , (6) 

where At is the time step. Analogously, basing on the equations for rotational motion 
(2), we define angular accelerations i?j and, therefore, two-fold time derivatives Aj(t ) 
(5), using principal torques Ki(to) and taking into account that f2j = AjA^~. So that the 
matrices A, at time to + At can be evaluated as follows 

A,(t + At) = Ai(to) + Mto)At + Mt )At 2 /2 + 0(At 3 ) . (7) 

And now we consider how to perform the second line 

s(t + At) = s(t ) + (s(t ) + s(t + At)) At/2 + 0(At 3 ) (8) 

of the velocity Verlet framework, where s denotes a spatial coordinate. There are no 
problems to pass this step in the case of translational motion, when s = fi and, therefore, 
for new translational velocities one obtains 

n(t + At) = n(t ) + (n(t ) + ri(t + At)) At/2 + 0(At 3 ) , (9) 

where fj(t + ^At) = ^-Fj(t + Z\t) and the forces Fi(t + /it) are calculated in the already 
defined new spatial configuration {rj(t + At), Aj(t + At)}. 

However, the difficulties immediately appear in the case of rotational dynamics, be- 
cause then second time derivatives s can depend explicitly not only on spatial coordinates 
s, associated with the rotational degrees of freedom, but also on generalized velocities s. 
For example, according to Euler equations (2), the principal angular accelerations depend 
on orientational variables via molecular torques and on angular velocities of molecules as 
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well. Then, choosing s = Aj, we obtain on the basis of the equations of motion (5) that 
Aj(t) = Aj(Aj(i), Aj(t)). In view of (8), this leads to a very complicated system of nine 
nonlinear equations per molecule with respect to the nine unknown elements of matrix 
Ai(t + At). It is worth to note that similar problems arise within the leapfrog and Beeman 
frameworks (see Appendix, where a rotational-motion version of the Beeman algorithm 
is derived). 

An alternative has been found in rotational- motion versions [21, 22] of the leapfrog 
algorithm. It has been assumed to associate the quantity s with the angular momentum 
li = A.fLi of the molecule in the laboratory system of coordinates, i.e., s = Zj, where 
Li = (Jif2{, J 2 /?2, J^l) = ifii and J is the diagonal matrix of principal moments of 
inertia. The rate of change in time of angular momentum is the torque, i.e., U = k t . 
Then equation (8) leads to a much more simple expression, 

h(t + At) = h(t ) + (ki(t ) + Hto + At)) At/2 + 0(At 3 ) , (10) 

and, therefore, new angular momenta are easily evaluated using the known torques k{ at 
times to and to + At. The corresponding values for principal angular velocities and first 
time derivatives of orientational matrices can be obtained, when they are needed, using the 
relations f2i(t +At) = Ai(t + At)li(t + At) and Ai(t +At) = Qi(t +At)Ai(t +At). 

Finally, we consider the third version of the velocity Verlet method. The idea consists 
in using angular velocities as independent parameters for describing the sate of the system 
in phase space. Then putting s = fii in (8) and taking into account Euler equations (2), 
we find the following result 

ak M 



2J, 

x 



K^t) + K l a (t + At) + (J? - J,) 

II 

(20^0)0^0) + nii^Afy + n^AQ}, + AfTpAW, 



The equations (11) constitute the system of maximum three nonlinear equations per 
molecule with respect to the same number of the unknowns AQ % a = ft l a {to + At) — f2 l a (to). 
The system (11) can be linearized, substituting initially AQ % a = in all quadratic terms, 
and solved in a quite efficient way by iteration. This is justified for At — > because then 
terms nonlinear in Af2 l a are small. 



From the mathematical point of view, all the three representations s = Aj, Zj or fii 
are completely equivalent, because the knowledge of an arbitrary quantity from the set 
(Aj, Zj, fl,j) allows us to determine the rest of two ones uniquely. In the case of numerical 
integration the pattern is different, because then the investigated quantities are evalu- 
ated approximately. The choice s = Aj can not be recommended for calculations due to 
its complexity. The case of s = Zj, corresponding to the angular- momentum version of 
the Verlet algorithm, is the most attractive in view of the avoidance of nonlinear equa- 
tions. Actual computations show, however, that the best numerical stability with respect 
to the total energy conservation is reached in the angular- velocity version (11) of the 
Verlet algorithm, when s = fi^. This fact can be explained taking into account that a 
kinetic part, \ Y^=x{J\^\ + + h&f ), of the total energy is calculated directly from 
principal angular velocities. At the same time, to evaluate angular velocities within the 
angular-momentum version the additional transformations J?j = J _1 AjZj with approxi- 
mately computed matrices Aj and angular momenta Zj are necessary. They contribute 
additional portions into the accumulated errors at calculations of the total energy. 

Shifting the initial time t to t + At, the integration process is repeated for a next 
time step. In such a way, step by step the dynamics of the system can be evaluated. 

B. Solving the rigidity problem 

Let us write an analytical solution for orientational matrices in the form 

P AfP 

Mt + At) = Y,A ( ?\t )^- + 0(At p+1 ) , (12) 
P =o P ] 

where A^\to) denotes the p-fold time derivative of Aj at time to- It can be shown easily 
from the structure of equation (3) that arbitrary-order time derivatives of the matrix 
constraint ®i(t) = Aj(t)A+(t) —1 = are equal to zero at a given moment of time, i.e., 
AjA+ + AjA+ = 0, AjA+ + 2AjA+ + AjA+ = and so on, when Aj is orthonormal. 
Therefore, if all the terms (P — > 00) of Taylor's expansion (12) are taken into account, that 
corresponds to the exact solution of equations of motion, and initially all the constraints 
are satisfied, 0j(t o ) = 0, they will be fulfilled at later times as well. 

In particular algorithms the expansion is truncated after a finite number of terms. 
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For example, the velocity- Verlet form (7) is restricted by quadratic terms (P = 2), in- 
volving truncation errors of order At 3 into the matrix elements of Aj. The same order 
of uncertainties will be accumulated in &i(t) at each time step, breaking the molecular 
structure, i.e., &i(to + At) = 0(At 3 ). In such a case the molecules collapse and can even 
be destroyed completely after a sufficient period of time. Therefore, the problem arises: 
how to modify the first line of the algorithm to achieve the exact rigidity for arbitrary 
times? 



1. Constraint-matrix scheme 

The usual way to reduce orientational matrices to orthonormal form lies in using 
the constraint technique. The main idea is simple. As far as the elements of orientational 
matrices are not independent, this requires, generally speaking, the necessity of introduc- 
ing additional forces which appear as a result of the constraints &i(t) = 0. These matrix 
constraints constitute, in fact, six independent scalar relations per molecule, namely, 



= a\l + a\ 2 2 + a\ 2 3 -1 = 0, 4>\ = a^a^ + a\ 2 a 22 + a\ 3 a 23 = , 



i l2 



l 22 



i 2 i 

aL - 1 



o 3 i + a m + a 



32 



°23 

;2 
*33 



0, 0*5 
1 = 0, <\>\ 



°ll a 31 ~l~ Q 19 a 3S> + Q 13 Q 



fl 2i a 31 



12"-32 



o , 



(13) 



t 23 u '33 



Then the corresponding constraint forces, acting on dynamical variables a l a/3 , are G l a/3 
— X)f=i ^\9<P\ I dd 1 ^ or in the matrix representation 



d = -AjAj 



A4 2A^ Aq 
\ A5 Aq 2 A3 



Aj , 



(14) 



where Aj are symmetric matrices of Lagrange multipliers. The matrices of constraint 
forces are now added in the equations of motion (5) and, as a consequence, the evaluation 
of matrix elements (7) is modified to 



A,(t + At) = Ai(io) + ki{t )At + Ai(t )At 2 /2 + Gi(t )At 2 /2 + 0{At 3 ) 



(15) 



In view of equations (14) and (15), to satisfy the conditions &i(t + At) = it is 
necessary to solve the system 0J(t o + At) = of six (/ = 1, . . . , 6) nonlinear equations 
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per molecule for six unknown Lagrange multipliers A; (to)- As usually, such a system 
is linearized and the unknowns are found by iteration. The iteration procedure can be 
initiated by putting A^ = in all nonlinear terms and the iterations converge rapidly at 
actual step sizes to the physical solutions A; (to) ~ At. The contributions of constraint 
forces into the matrix evaluation (15) are of order At 3 , i.e., the same order as truncation 
errors of the basic algorithm (7), but the rigidity is now fulfilled perfectly for arbitrary 
times in future. It is worth to remark that the constraint forces introduced should be 
treated as pseudo forces, because they depend on details of the numerical integration in a 
characteristic way and disappear if the equations of motion are solved exactly, i.e., when 
At -> 0. 

2. Rotational-matrix scheme 

Fortunately, the cumbersome procedure of solving nonlinear equations to preserve 
the molecular rigidity can be avoided in our approach using the fact that actual algorithms 
are accurate to a finite order only in time step. In view of equalities (3) and (5), the 
evaluation (7) can be presented in a more compact form, 

Ai(i + At) = Di(t , At)Ai(t ) + 0(At 3 ) , (16) 

where 

T>i(t , At) —1 + QMAt + (pitto) + fi?(t )) At 2 /2 (17) 

are evolution matrices. Let the rigidity has been satisfied at time to, i.e., 0j(io) — 
A i (t )A+(t ) -1 = 0. Then 0,(t o + At) = T>i(t , At)T>? (t , At) - I = 0(At 3 ) or, in 
other words, the matrices Dj are not orthonormal. 

The simplest way to present the evolution matrices and, as a consequence, the ori- 
entational matrices in orthonormal form lies in the following. Taking into account that 

= W(J?j) — fii I, where Qi = \J Q\ + Q\ + f2% is the magnitude of the angular 
velocity and W(i?j) is a symmetric matrix with the elements fi l a &p, we rewrite (17) as 

T>i(t , At) = (1 - ^ 2 (t )Z\t 2 /2)I + ^i{Oi{to))At 2 /2 + Ui(t )At , (18) 
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where fii(to) is an antisymmetric matrix of type (4), constructed on the mean value 
f2i(t ) = Oi(t ) + f2i(t )At/2 of the angular velocity for the i-th molecule during the 
time interval [to, to + A]. It is easy to see that replacing i?j by i?j in (18), we introduce 
the error of order At 3 . Moreover, taking into account that 

At = — ^ — '- + 0{At 3 ) , — = 1 ^ — i + C(Z\t 4 ) , (19) 

we adjust (18) to the form 

V t (t , At) = IccxfaAt) + 1 ~ <X _^ 2iAt) W(n i ) + ^^ U, = eMVMAt) . 

(20) 

Let us expand the matrix X>j(t , At) into the Taylor's series with respect to At. Then 
it can be verified easily that each elements of this matrix coincides with the corresponding 
element of T)i(to,At) (17) up to the second order in At inclusively. Higher order terms, 
being associated with time derivatives of angular accelerations, are not taken into account 
within the velocity Verlet framework and they can merely be omitted without loss of the 
precision. Therefore, the matrices / D i (to,At) (20) and T>i(t ,At) (17) differ between 
themselves by terms of order At 3 or higher that is completely in the self-consistency with 
truncation errors of the algorithm considered. However, the main advantage of using T> \, 
instead of Dj, in the evaluation 

A,(t + At) = T>i(t , At) AM + 0(At 3 ) , (21) 

of orientational variables consists in the fact that the matrix 'D i (to,At) is orthonormal, 
\.e.,T> i (to,At)T>f{to,At) = I and then @i(t + At) = T>i(t , At) g D?(t , At) -I = 0. As it 
follows from the structure of eq. (20), the matrix X>j(t , At) defines the three-dimensional 
rotation on angle Q^At around the axis directed along vector In such a way, the rigid 
structures of molecules can be reproduced exactly at each time step of the integration. 

IV. NUMERICAL TESTS AND DISCUSSION 

We now test our matrix method on the basis of simulations on a TIP4P model 
[24] of water. This method was used by us previously [25] investigating a Stockmayer 
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fluid of point dipoles. In the TIP4P model the water molecule consists of four sites, 
M — 4. We have used a system of N = 256 molecules and the interaction site reaction 
field geometry [26]. Intersite components of the TIP4P potential were cut off and shifted 
to zero at point of the truncation to avoid the system energy drift associated with the 
passage of the sites through the surface of the cut-off sphere. The cut-off radius was 
half the basic cell length. The MD simulations were performed in the microcanonical 
ensemble at a density of 1 g/cm 3 and at a temperature of 298 K. The numerical stability 
of solutions to the equations of motion was identified in terms of relative fluctuations, 
E(t) = \J ' {{E — (E) t ) 2 )t/ (E)t, of the total energy E of the system during time t. 

We have made a comparative test carrying out explicit MD runs using our angular- 
velocity Verlet integrator (eq. (11)) within constraint- and rotational-matrix schemes 
(eqs. (15) and (21), respectively), as well as the implicit quaternion leapfrog algorithm 
[22], the pseudo-particle formalism [23] and the atomic-constraint technique [13]. The 
runs were started from an identical well equilibrated configuration. All the algorithms 
required almost the same computer time per step (96% being spent to evaluate pair in- 
teractions). For the purpose of comparison the quaternion integration with the Gear 
predictor-corrector algorithm of fifth order [15, 16] has been considered as well. At least 
two corrector steps were used to provide an optimal performance of the predictor-corrector 
scheme and, as a consequence, twice or more larger computer time was taken in this case 
than that is normally necessary. 

The results obtained for relative total energy fluctuations as functions of the length of 
the simulations at four fixed step sizes, At = 1, 2, 3 and 4 fs, are presented in fig. 1 (water 
is most commonly simulated with a step size of order 2 fs [27]). At small time steps, 
At < 1 fs, all the approaches exhibited similar equivalence in the energy conservation 
(subset (a) of fig. 1), except the Gear algorithm which produced much more accurate 
trajectories. But the Gear algorithm begins to be unstable already at At — 1 fs and leads 
to the worst results for At > 1.5 fs (see, as an example, the case At — 2 fs, subset (b)). 
Somewhat better stability is observed in the leapfrog and pseudo-particle approaches. 
However, at moderate and great time steps, At > 2 fs (figs. 1 (b)-(d)), the results are 
rather poor, especially in the case of the leapfrog scheme. The best numerical stability 
has been achieved with the atomic-constraint algorithm and our matrix method, which 
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conserve the total energy approximately with the same accuracy up to At — 3 fs. It can 
be seen easily that the matrix method works better within the rotational-matrix scheme, 
so that there is no need to use the complicated constraint-matrix procedure. Quite a few 
iterations (the mean number of iterations per molecule varied from 3 to 5 at At = l-j-4 fs) 
was sufficient to obtain solutions to the system of quadratic equations (11) with a relative 
iteration precision of 10~ 12 . This contributes a negligible small portion additionally into 
the total computer time. 

To demonstrate that the exact reproduction of molecular rigidity is so important, we 
have also integrated the equations of motion in a situation (eq. (7)) when no additional 
normalization and orthogonalization of principal-axis vectors are used. In this case the 
total energy fluctuations increased drastically with increasing the length of the runs at 
arbitrary time steps (see, for instance, the corresponding curve in fig. 1 (a)). The same 
words can be said in the case when no quaternion renormalization is applied along the 
leapfrog trajectories. This is so because in the free-normalization regime, the structure of 
molecules is broken that leads to an unpredictable discrepancy in the calculation of po- 
tential forces and significant deviations of the total energy. We have also established that 
the numerical stability is very sensitive to the way of how the quaternion renormalization 
is performed. In particular, the energy conservation can be somewhat improved if the 
quaternions are renormalized inside the iterative loop of the implicit leapfrog integrator 
rather than at the end of each time step only, as was originally proposed [22]. 

No shift of the total energy has been observed for the atomic-constraint and matrix 
approaches over a length of 10 ps at At < 3 fs. Instead, it oscillates around a stable 
value of E = —33.6 kJ/mol. To reproduce the features of microcanonical ensembles 
quantitatively, it is necessary for the ratio r = S/T of total energy fluctuations to fluc- 
tuations T of the potential energy to be no more than a few per cent. For the system 
under consideration T 0.56%, so that, for example, the level £ = 0.03% will corre- 
spond to T ~ 5% that is still acceptable for precision calculations. The ratios T, obtained 
within various approaches at the end of 10 ps runs, are plotted in fig. 2 as dependent 
on the time increment. The results of fig. 2 show that a level of T = 5% is achieved at 
the time steps 1.2, 1.4, 3.0 and 4.0 fs within the leapfrog, pseudo-particle, matrix and 
atomic approaches, respectively. Therefore, the last two methods allow a step size more 
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than twice larger than the pseudo-particle and leapfrog algorithms. The functions T(At) 
can be interpolated with a great accuracy as C At 2 + C'At 3 with the coefficients C w 
0.28 and 0.30 % fs~ 2 , C 0.01 and 0.10 % fs~ 3 for the atomic and matrix approaches, 
respectively. The characteristic square growth of V at small time steps is completely in 
line with 0(At 2 ) order of global errors for the algorithm considered. 

It is worth to underline that analyzing the system over a significantly shorter time 
period, of order 1 ps say (as was done by Ahlrichs and Brode [23]), one may come to a very 
misleading idea about the energy conservation. We can see clearly from fig. 1 that such a 
simulation period (corresponding to 1000, 500, 333 and 250 time steps at At = 1, 2, 3 and 4 
fs, respectively) is quite insufficient to give a realistic pattern on global errors accumulated 
in the total energy. And only beginning from lengths of order 10 ps, we are entitled to 
formulate true conclusions on the numerical stability. These lengths are sufficiently long 
to observe an appreciable modification of the system. For instance, during 10 ps even long- 
lived dipole moment correlations vanish completely [28]. Moreover, the phase trajectories 
of 10 ps long are also sufficient, as a rule, to reproduce thermodynamic, structure and 
other properties of water with a reliable statistical accuracy. The investigation of some 
collective effects, such as dielectric relaxation, may require extremely long simulations (up 
to 1000 ps [28]) to reduce statistical noise. As a result, even the best algorithms may not 
provide a required numerical stability. In such a situation, we can merely slightly rescale 
the velocities of particles when the total energy has exceeded an allowed level. Obviously, 
the investigated quantities will be little affected by this rescaling if it is applied not more 
frequently than after a period of time during which the correlations have significantly 
decayed. 

In view of the results obtained in this section, we can conclude that the method 
proposed appears to be the most efficient among all known algorithms deriving within 
the molecular framework and can be considered as a good alternative to the cumbersome 
atomic technique. The fact that our molecular Verlet algorithm conserves the total energy 
at great step sizes somewhat worse than the atomic Verlet algorithm results from the 
introduction of velocities. As far as velocities appear explicitly, the angular accelerations 
begin to be velocity dependent. Further, the angular velocities are calculated with one step 
errors of order 0(At 3 ) and the same order of uncertainties will be presented simultaneously 
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in angular accelerations. This, in its turn (see eq. (11)), leads to additional terms of orders 
0(At 4 ) and 0(At 3 ) in the truncation and global errors, respectively, for angular velocities 
and, as a consequence, for the total energy. That is why in the case of rotational motion the 
coefficient C corresponding to the velocity Verlet differs significantly from that obtained 
for the usual (free of velocities) Verlet algorithm. At the same time, the corresponding 
values of C are practically equal between themselves, and, therefore, we may stay about 
the equivalence of the both algorithms with respect to the main term of global errors. 

The pointed out above minor disadvantage is compensated, however, by a much more 
major advantage of our method with respect to the atomic scheme in that the velocity 
Verlet algorithm allows to perform simulations in canonical ensembles. As is well known 
[22], thermostat calculations can be carried out with significantly greater step sizes than 
those used in microcanonical ensembles. A thermostat version of the velocity Verlet 
algorithm for rotational motion will be studied in a separate investigation. 

V. CONCLUSION 

We have shown that the difficulties in numerical integration of rigid polyatomics can 
be overcame using an alternative approach. In our singularity free scheme, orientational 
matrices were used to represent the rotational degrees of freedom of the system. Although 
this introduces extra equations per molecule and the lack of independence for the matrix 
elements, but presents no numerical difficulties. An elegant procedure, built directly into 
the Verlet algorithm, has allowed to perfectly fulfil the rigidity of molecules at each step 
of the trajectory without any additional efforts and loss of precision. Avoidance of the 
necessity to solve complex nonlinear equations for preservation of the molecular rigidity 
should be a benefit of the matrix method with respect to the atomic-constraint approach. 

We have demonstrated on the basis of actual calculations that the matrix method leads 
to results comparable in efficiency with the cumbersome atomic-constraint technique. The 
advantages of the matrix scheme are that it can be implemented for arbitrary rigid bodies, 
extended to a thermostat version and realized in MD programmes in a more simple way. 
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Appendix 



We now consider the question of how to adopt our matrix scheme to integrate the equations 
of motion within the Beeman framework. According to the usual Beeman algorithm [20], the 
translational positions and velocities of molecules are evaluated as 

n(to + At) = n(t ) + ri(to)At + [§fi(*o) - b^(*o - At)]At 2 + 0{At A ) , 

(Al) 

ri(t + At) = fi(t ) + [ifi(*o + At) + |fi(*o) - \ri(t Q - At)]At + 0{At 3 ) . 

The order of truncation errors in coordinates increases to four because the expression [§ri(io) — 
ifiOo - At)]At 2 can be reduced to the form fi(t )At 2 /2 + ir^t^Afi/Q + 0(At A ) with the 
estimation r'j(io) = [f""i(io) — f""i(*o — At)]/ At + O(At) of superaccelerations. The fractions in 
the second line of eq. (Al) are obtained in such a way to provide the third order of truncation 
errors in velocities and to satisfy exactly the Stormer central difference approximation [16, 29] 
of accelerations 

s(t + At) = -s(t - At) + 2s(t ) + s(t )At 2 + <D(At A ) (A2) 

with s = rj. Acting in the spirit of the Beeman framework, we can write analogous to (Al) 
equations for orientational matrices and angular velocities. The result is 



Ai(t + Aty^i(t ) + Ai{t )At + [§Ai(t ) - iMto ~ At)]At 2 

-ilM^Mito) - iA^o - At)A t (t - At)]At 2 + 0(At A ) , 



(A3) 



&J n+1 \t + At) = fl^to) + ^ 



lK(t + At) + lK(t) - lK(t - At) + (jp - j 7 ) 

(|4 (n) (t + ^)«* (n) (t + At) + §l2j(t )^(to) - g^fo - ^t)^(*o - 4t 



(A4) 



where the symmetric constraint matrices A, (to) ~ At 2 are found from the constraint relations 
Aj(to + At)Af(to + = I, whereas new values ^(to + ^) for principal components of the 
angular velocities can be computed by iteration (n = 0, 1, . . .) taking J2^°^(£o + At) = ^ l a (to) 
as initial guesses. 

A rotational-matrix scheme can be derived within the Beeman method as follows. Consider 
first a more general procedure for the orthonormalization of orientational matrices, which will 
be valid for integrators of arbitrary order in truncation errors. Let the algorithm applied uses 
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Taylor's expansion (12) for the time evaluation of orientational matrices. Then the evolution 
matrices can be cast as 

B l (t ,At)=I + J2^ ) (to)^- , (A5) 
P =i P- 

where d| p) = a\ p } A+, or more explicitly: dJ 1} = n i5 Df ) = + O?, Df ) = ft; + 20^ + 
ftifti + , Df } = ft; + 3ft;ft; + ft;ft; + 3ft;ft; + 3ft;ft 2 + 2ft;ft;ft; + ft 2 ft; + ftf and so on. 
A rotational-matrix counterpart of (A5) we find in the orthonormal form 

Vi(t ,At) =exp(E;=iHl p) (t )^) , (A6) 

where are unknown antisymmetric matrices, i.e., H- p ^ = — Hf \ and expand the exponent 
(A6) into the Taylor series at At — > 0. It is obvious that -D;(io> At) and T>i(to, At) will be iden- 
tical at P — > oo, if all their matrix coefficients, corresponding to the same powers p = 1, 2, . . . , P 
of Z\t, are equal between themselves. This condition leads to a recursive procedure with the 
solutions = ft;, Hj 2) = fli, Hf } = ft; + ±(ft;ft; - ft;ft;), hJ 4) = 4 + - 

and so on. The Beeman approach is accurate to third order in coordinates (P = 3), i.e., 
[§A;Oo) - gA»(to - At)]At 2 = Ai(t )At 2 /2 + Ai(t )At 3 /6 + C(Z\t 4 ), where the super accelera- 
tions Aj(t ) = [Aj(io) — Aj(t — Z\t)]/Z\i + C(Z\t). Similarly we can estimate angular superac- 
celerations, ft;(to) = [ftj(i ) — ftj(io)(£o ~~ At)]/ At, and obtain in this case 

Ttf{h,At) = exv(n i (t Q )At+[\tuW . 

(A7) 

Putting P = 2 in eq. (A6) yields the result Z>i(£ , At) = exp(fti(t )At + ±ft;(i )Z^ 2 ). As was 
expected, this is completely in line with the result (20) performed in Sec. 3 for the velocity Verlet 
algorithm on the basis of intuitive grounds. 

It is worth mentioning that approximation (A2) is used directly for evaluation of spatial 
coordinates in the usual Verlet algorithm [17, 23]. As can be verified, any trajectory produced 
by the velocity Verlet algorithm satisfies equation (A2) at s = {r;,A;} even if constraint- 
or rotational-matrix schemes are used. The fact that the trajectory s(t) can be generated 
with the same fourth-order local errors by lower-order equations (6) and (7) (or (21)) results 
from a fortunate cancellation of truncation errors arising in coordinates and velocities during 
two neighbour time steps. Note, however, that the usual Verlet algorithm, its velocity version 
and Beeman method are not equivalent, because they differ between themselves by the main 
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term of fourth-order uncertainties in coordinates and calculate one-step velocities in a different 
manner. For example, evaluating velocities within the Beeman algorithm, it is assumed that the 
accelerations are slow variables on time scales of 2 At. If this criterion is not satisfied, the main 
term 0(At s ) of truncation uncertainties in velocities and, as a result, the main coefficient C in 
global errors for the total energy may increase in a characteristic way. This prediction has been 
confirmed by our computer simulations on the TIP4P water. Therefore, the Beeman algorithm 
can be applied for systems with sufficiently smooth interparticle potentials only. 
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Figure captions 



Fig. 1. The relative total energy fluctuations as functions of the length of the simu- 
lations on the TIP4P water, obtained within various techniques at four fixed time steps, 
namely, 1 fs (a), 2 fs (b), 3 fs (c) and 4 fs (d). 

Fig. 2. The ratios of the total energy and potential energy fluctuations as dependent 
on the step size, observed for various approaches in the simulations of the TIP4P water 
at the end of 10 ps runs. 
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